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Abstract 

A discrete model describing defects in crystal lattices and having the standard linear anisotropic 
elasticity as its continuum limit is proposed. The main ingredients entering the model are the 
elastic stiffness constants of the material and a dimensionless periodic function that restores the 
translation invariance of the crystal and influences the Peierls stress. Explicit expressions are 
given for crystals with cubic symmetry: sc, fee and bcc. Numerical simulations of this model 
with conservative or damped dynamics illustrate static and moving edge and screw dislocations 
and describe their cores and profiles. Dislocation loops and dipoles are also numerically observed. 
Cracks can be created and propagated by applying a sufficient load to a dipole formed by two edge 
dislocations. 

PACS numbers: 61.72.Bb, 5.45.-a, 82.40.Bj, 45.05. +x 
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I. INTRODUCTION 



The advances of electronic microscopy allow imaging of atoms and can therefore be used 
to visualize the core of dislocations [Jj, |2| , cracks 3} and other defects that control crystal 
growth and the mechanical, optical and electronic properties of the resulting materials 0]. 
Emerging behavior due to motion and interaction of defects might explain common but 
poorly understood phenomena such as friction nil. Defects can be created in a controlled 

n 

way by ion bombardment on reconstructed surfaces p , which allows the study of effectively 
two dimensional (2D) single dislocations and dislocation dipoles. These dislocations are 
effectively 2D because the surface 'floats' on the 3D crystal 7L Other defects that are 
very important in multilayer growth are misfit dislocations jg, 2, At the nanoscale, 
many processes (for example, dislocation emission around nanoindentations involve 
the interaction of a few defects so close to each other that their core structure plays a 
fundamental role. To understand them, the traditional method of using information about 
the far field of the defects (extracted from linear elasticity) to infer properties of far apart 
defects reaches its limits. The alternative method of ab initio simulations is very costly 
and not very practical at the present time. Thus, it would be interesting to have systematic 
models of defect motion in crystals that can be solved cheaply, are compatible with elasticity 
and yield useful information about the defect cores and their mobility. 

To see what these models of defects might be like, it is convenient to recall a few facts 
about dislocations. Consider for example an edge dislocation in a simple cubic (sc) lattice 
with a Burgers vector equal to one interionic distance in gliding motion, as in Fig. 
The atoms above the xz plane glide over those below. Let us label the atoms by their 
position before the dislocation moves beyond the origin. Consider the atoms (xo, —a/2,0) 
and (xo,a/2,0) which are nearest neighbors before the dislocation passes them. After the 
passage of the dislocation, the nearest neighbor atoms are (x , —a/2, 0) and (x — a, a/2,0). 
This large excursion is incompatible with the main assumption under which linear elasticity 
is derived for a crystal structure the deviations of ions in a crystal lattice from their 
equilibrium positions are small (compared to the interionic distance), and therefore the ionic 
potentials entering the total potential energy of the crystal are approximately harmonic. One 
obvious way to describe dislocation motion is to simulate the atomic motion with the full 
ionic potentials. This description is possibly too costly. In fact, we know that the atomic 
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FIG. 1: (Color online) Deformed cubic lattice in the presence of an edge dislocation for the piecewise 
linear g(x) of Eq. @ with a = 0.24. 



displacements are small far from the dislocation core and that linear elasticity holds there. 
Is there an intermediate description that allows dislocation motion in a crystal structure and 
is compatible with a far field described by the corresponding anisotropic linear elasticity? 

If we try to harmonize the continuum description of dislocations according to elasticity 
with an discrete description which is simply elasticity with finite differences instead of dif- 
ferentials, we face a second difficulty. The displacement vector of a static edge dislocation 
is multivalued. For example, its first component is u\ = a(27r) _1 [tan~ 1 (?//x) + xy/(2(l — 
v){x 2 + y 2 ))} for the previously described edge dislocation [v is the Poisson ratio) I n 
elasticity, this fact does not cause any trouble because the physically relevant strain tensor 
contains only derivatives of the displacement vector. These derivatives are continuous even 
across the positive x axis, where the displacement vector has a jump discontinuity [ui] = a. 
If we consider a discrete model, and use differences instead of differentials, the difference of 
the displacement vector may still have a jump discontinuity across the positive x axis. 

The previous difficulties have been solved in a simple discrete model of edge dislocations 
and crowdions called the IAC model (interacting atomic chains model) proposed and studied 
by A.I. Landau and collaborators 1311 . A similar model for screw dislocations in bcc crystals 
was proposed earlier by H. Suzuki [14]. In the equations for the IAC model, the differences 
of the displacement vector are replaced by their sines. Unlike the finite differences, these 
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sine functions are continuous across the positive x axis. Moreover, the equations remain 
unchanged if a horizontal chain of atoms slides an integer number of lattice periods a over 
another chain. Taking advantage of its simplicity, we have recently analyzed pinning and 
motion of edge dislocations in the IAC model |15| . 

In this paper, we propose a top-down approach to discrete models of dislocations in cubic 
crystals. Let us start with a simple cubic lattice having a unit cell of side length a. Firstly, we 
discretize space along the primitive vectors defining the unit cell of the crystal: x = x\ = la, 
y = x 2 — ™, z = £3 = na, where /, m and n are integer numbers. We shall measure the 
displacement vector in units of a, so that Ui(x,y,z,t) = aui(l,m,n;t) and Ui(l,m,n;t) is 
a nondimensional vector. Let D~j~ and D~ represent the standard forward and backward 
difference operators, so that Dfu{(l, m, n; t) — ± ± 1, m, n; t) — m, n; t)}, and so on. 
We shall define the discrete distortion tensor as 
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(?) 



giD+m), (l) 



where g(x) is a periodic function of period one satisfying g(x) ~ x as x — > 0. In this paper, 
we shall use the odd continuous piecewise linear function: 

ftC ^ I iC I 2 5 
(1-2^(1-2.) l_ a<x< l ^ 

4a ' 2 2 ' 

which is periodically extended outside the interval (—1/2, 1/2) for a given a, < a < 1/2. 
Note that g(x) is symmetrical if a = 1/4 and that the interval of x in which g'(x) < 
widens with respect to that in which g'(x) > as a increases. Numerical simulations of 
the governing equations for a 2D edge dislocation show that the Peierls stress decreases as 
a increases; see Fig. which will be further commented later on. This means that the 
dislocation is harder to move if a decreases, i.e., if the interval of x in which g'(x) < 
shrinks with respect to that in which g'(x) > 0. The parameter a can be selected so as to 
agree with the observed or calculated Peierls stress of a given crystal. 
Secondly, we replace the strain tensor in the strain energy by 



e, = \ («,?> + wf) = 9&M+jVtej. . (3) 
Summing the strain energy over all lattice sites, we obtain the potential energy of the crystal: 



V({ui}) = a 3 W(l,m,n;t). (4) 

l,m,n 
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FIG. 2: (Color online) Peierls stress in dimensionless units for a 2D edge dislocation in a sc crystal 
with the stiffnesses of tungsten as a function of the parameter a in the periodic function g(x). 

W(l,m,n;t) = W({ui}) = \^i^%fikh (5) 
Cijki = A SijSia + /i (S ik 5ji + 8u8 jk ) + 2(C 44 - p) 

SiiSij5i k 5u - 52iS 2 jS 2 kS2i - S3i8 3j S 3k 8 3 i , (6) 



8ik8ji + 8u5jk 



in which summation over repeated indices is understood. Here, A = Ci 2 , /i = (Cn — Ci 2 )/2, 
where Cij are the stiffness constants of a cubic crystal. If C44 = fi, the strain energy is 
isotropic and A and fi are the usual Lame coefficients. 

Next, we find the equations of motion. In the absence of dissipation and fluctuation 
effects, they are 



4 •■ /, A 1 dV({u k }) ( 

pa Ui(l,m,n;t) = — ^— -, 7 

a 0Ui(L, m, n; t) 



or, equivalently, 



MuS,™,n;t) = - JL— - £ W(l',m',n';t), (8) 
M = pa 2 . (9) 

Here iij = d 2 Ui/dt 2 , M has units of mass per unit length (because p is the mass density) and 
the displacement vector is dimensionless, so that both sides of Eq. (jHj) have units of force 
per unit area. We show in Appendix |X] that Eq. (jHJ) is equivalent to the following spatially 
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discrete equations 



Mui 



52D-[c m g'(Df Ui )g(D+u k )}. 



(10) 



j,k,l 



To nondimensionalize these equations, we could adopt C44 as the typical scale of stress and 



t = JM/C44 = ayp/C*44 as the time scale. The resulting equations are the same ones with 
M = 1 and c^/C^ instead of Cijki- Let us now restore dimensional units to this equation, 
so that iti(x,y,z) = aUi(x/a,y/a, z/a), then let a — > 0, use Eq. © and that g(x) ~ x as 



Thus the discrete model with conservative dynamics yields the Cauchy equations for elastic 
constants with cubic symmetry provided the components of the distortion tensor are very 
small (which holds in the dislocation far field). Equations of motion with dissipation and 
fluctuation terms can be obtained by writing a quadratic dissipative function which, in the 
isotropic case, yields the usual fluid viscosity terms and using the fluctuation-dissipation 
theorem \)\ . 

In the rest of the paper, we describe dislocation motion in simple cubic crystals and 
extend our discrete elastic equations to the case of fee and bec lattices. 

II. DISLOCATION MOTION IN SC CRYSTALS 

In this Section, we shall find numerically pure screw and edge dislocations of our dis- 
crete model for sc symmetry and discuss their motion under appropriate applied stresses. 
In all cases, the procedure to obtain numerically the dislocation from the discrete model 
equations is the same. We first solve the stationary equations of continuum elasticity with 
appropriate singular source terms to obtain the dimensional displacement vector u(x, y, z) = 
(ui(x, y, z), m 2 (x, y, z),u 3 (x, y, z)) of the static dislocation under zero applied stress. This dis- 
placement vector yields the far field of the corresponding dislocation for the discrete model, 
which is the nondimensional displacement vector: XJ(l,m, n) = u(la,ma,na)/a. We use the 
nondimensional static displacement vector U(/,m, n) in the boundary and initial conditions 
for the discrete equations of motion of the discrete model. Later in the Section, we shall 
show numerical results corresponding to the interaction of edge dislocations and the opening 
of a crack. 




x — > 0. Then we obtain the equations of linear elasticity [16], 
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FIG. 3: (Color online) Screw dislocation for the piecewise linear g{x) of Eq. @ with a = 0.24. 



FIG. 4: (Color online) Superimposed zooms of the moving core of a screw dislocation. 
A. Screw dislocations 

The continuum displacement field of a dislocation, u = (€t 1? u 2 , U3), can be calculated as 
a stationary solution of the anisotropic Navier equations with a singularity oc r _1 at the 
dislocation core and such that J c (c?x- V)u = — b, where b is the Burgers vector and C is any 
closed curve encircling the dislocation line A pure screw dislocation along the z axis with 
Burgers vector b = (0, 0, b) has a displacement vector u = (0, 0, 1X3(2;, y)) I}. Then the strain 
energy density © becomes W = C44 |Vm3| 2 /2, and the stationary equation of motion is 
A£t 3 = 0. Its solution corresponding to a screw dislocation is u^{x, y) —b (27r) _1 tan _1 (|//a;) 
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[lj]. The same symmetry considerations for Eq. (|10j) yield the following discrete equation for 
the z component of the nondimensional displacement u 3 (l,m;t): 

Mu 3 = C u {D^[g(Dtu 3 )g'(Dtu 3 )]+D^[g(D+u 3 )g'(D+u 3 )}}. (12) 

Numerical solutions of Eq. (JT2*j) show that a static screw dislocation moves if an applied 
shear stress surpasses the static Peierls stress, \F\ > F cs , but that a moving dislocations 
continues do,n g so until the applied shear stress falls below a lower threshold F ci (dynamic 
Peierls stress); see Ref. |15j] for a similar situation for edge dislocations. To find the static 
solution of this equation corresponding to a screw dislocation, we could minimize an energy 
functional. However, it is more efficient to solve the following overdamped equation: 

(3u 3 = C M {D^[g(D+u 3 )g'(D+u 3 )}+D^[g(D+u 3 )g'(D+u 3 )}}. (13) 

The stationary solutions of Eqs. ()12j) and (fT3*j) are the same, but the solutions of (fT3*|) relax 
rapidly to the stationary solutions if we choose appropriately the damping coefficient (3. We 
solve Eq. |T3j) with initial condition u 3 (l,m;0) = U 3 (l,m) = b (2ira)^ 1 tan _1 (m//), and with 
boundary conditions u 3 (l, m; t) = U 3 (l, m) + Fmat the upper and lower boundaries of our 
lattice. At the lateral boundaries, we use zero-flux Neumann boundary conditions. Here F 
is an applied dimensionless stress with \F\ < F cs ; to obtain the dimensional stress we should 
multiply F by C44. For such small stress, the solution of Eq. (|13p relaxes to a static screw 
dislocation u 3 (l, m) with the desired far field. If F — 0, Figure El shows the helical structure 
adopted by the deformed lattice (l,m,n + u 3 (l,m)) for the piecewise linear g(x) of Eq. (J2J) 
with a = 0.24. The numerical solution shows that moving a dislocation requires that we 
should have g'(D+u 3 ) < (with either j — 1 or 2) in (|12|) or (jT5j) at its core [15]. This is 
harder to achieve as a decreases. 

The motion of a pure screw dislocation is somewhat special because its Burgers vector 
is parallel to the dislocation line. Any plane containing the Burgers vector can be a glide 
plane. Under a shear stress F > F cs directed along the y direction, a screw dislocation moves 
on the glide plane xz. A moving screw dislocation has the structure of a discrete traveling 
wave in the direction x, with far field u 3 (l — ct, m) + Fm; c = c(F) is the dislocation speed; 



see Figure EJ T 
as discussed in 



lis is similar to the case of edge dislocations in the simple IAC model 



[13], 



15j where the details of the analysis can be looked up. 
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B. Edge dislocations 



To analyze edge dislocations in the simplest case, we consider a isotropic cubic 
crystal (C 44 = (Cn — Ci 2 )/2) with planar discrete symmetry, so that u(l,m;t) = 
(ui(l, m; t), u 2 (l, m; t), 0) is independent of z = na. 

To find the stationary edge dislocation of the discrete model, we first have to write the 
corresponding stationary edge dislocation of isotropic elasticity. An edge dislocation directed 
along the z axis (dislocation line), and having Burgers vector (6,0,0) has a displacement 
vector u = (u x (x, y), u 2 (x, y), 0) with a singularity oc r _1 at the core and satisfying J c (<ix ■ 
V)u = —(b, 0, 0), for any closed curve C encircling the z axis. It satisfies the planar stationary 
Navier equations with a singular source term: 

1 



Au + 



l-2u 



V(V ■ u) 



-(0,6,0)5(r). 



Here r = y/x 2 + y 2 and v = A/[2(A + fi)] is the Poisson ratio; cf. page 114 of Ref. jl^ . 
appropriate solution is 



(14) 



The 



Ul 



b_ 

2^ 



tan -1 [ - I + 



Uo = 

2tt 



1 - 2u 



xy 



2(1 - u)(x 2 + y 2 ) 



In 



x 



+ 



(15) 



4(l-i/) \ b 2 J 2(1 - v)(x 2 + y 2 ) 
cf. Ref. |J, pag. 57. 

Eqs. (115)) yield the nondimensional static displacement vector U(/,m) = 
(ui(la,ma)/a,U2(la,rna)/a,0), which will be used to find the stationary edge dislocation 
of the discrete equations of motion. For this planar configuration, the conservative equa- 
tions of motion (fTUj) become 



Mu x = C xl Dl\ g (Dtu x )tf(D+u x )\ + C V iDi\ g (Dtvv)g\D$ui)\ 
+C 44 D 2 {[g(D+u x ) + g(D+u 2 )\ g'(D+u x )}, 

Mu 2 = C xl D 2 [g(D+u 2 )g'(D+u 2 )} + C X2 D 2 [g(D+u x )g'(D+u 2 )} 
+C U D^{[g{D+u 2 ) + g(D+u x )} g'(D+u 2 )}. 



(16) 
(17) 



To find the stationary edge dislocation corresponding to these equations, we set C 44 = 
(Cn — C X2 )/2 (isotropic case), and replace the inertial terms Mu x and Mu 2 by (3ii x and f3ii 2) 
respectively. The resulting overdamped equations, 



Pin = C n Dt\g(Dtu 1 )g , (D+u l )]+CuDt\g(D+u 2 )g , (,Dtu 1 )] 



FIG. 5: (Color online) Edge dislocation for the piecewise linear g(x) of Eq. © with (a) a = 0.27, 
(b) a = 0.29, (c) a = 0.32. 

+CuD 2 {[g(D+ Ul ) +g(Dfu 2 )]g\D+u 1 )}, (18) 
(3ii 2 = CnDMD+v*) g'(D+u 2 )} + C 12 D 2 lg(D+ Ul ) g'(D+u 2 )] 

+C u D^{\g(D+u 2 )+g(D+ Ul ))g'(Dtu 2 )}, (19) 

have the same stationary solutions as Eqs. (jTfij) and (JT7J). We solve Eqs. (fTHJl and (fEJJ) 
with initial condition u(Z, to; 0) = U(Z, to) given by Eqs. (|15|h and with boundary conditions 
u(Z, to; t) = U(Z, to) + (Fm, 0, 0) at the upper and lower boundaries of the lattice (F is a di- 
mensionless applied shear stress; recall that the displacement vector in the discrete equations 
is always dimensionless) . If \F\ < F cs (F cs is the static Peierls stress for edge dislocations), 
the solution of Eqs. ([THjl and (fTTIjl relaxes to a static edge dislocation (ui(l, to), u 2 (l, to), 0) 
with the appropriate continuum far field. 

In our numerical calculations of the static edge dislocation, we use the elastic constants of 
tungsten (which is an isotropic bcc crystal), Cn = 521 GPa, C\ 2 = 201 GPa, C44 = 160 GPa 
(Cn = C12 + 2C44) [2|. This yields v = 0.278. Figure El shows the structure adopted by the 
deformed lattice (/ + ui(l,m),m + U2{l,m)) when v = 0.278. for the asymmetric piecewise 
linear function g(x) with three different values of a. The profiles of the displacement vector 
are shown in Figures |H1 and [7| The glide motion of edge dislocations occurs on the glide 
plane defined by their Burgers vector and the dislocation line, and in the direction of the 
Burgers vector. In our case, a shear stress in the direction y, will move the dislocation in the 
direction x. For conservative or damped dynamics, the applied shear stress has to surpass 
the static Peierls stress to depin a static dislocation, and a moving dislocation propagates 
provided the applied stress is larger than the dynamic Peierls stress (smaller than the static 
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FIG. 6: (Color online) Profile of U\(l,m) for an edge dislocation with the piecewise linear g(x) of 




FIG. 7: (Color online) Profile of 112(1, m) for an edge dislocation with the piecewise linear g(x) of 
Eq. © when (a) a = 0.24, (b) a = 0.27, (c) a = 0.29, (d) a = 0.32. 

one) Q]. Numerical solutions show that the dimensionless Peierls stress depends on a in 
Eq. (j2J) as shown in Fig. [21 As the interval of x for which g'(x) < shrinks (which occurs as 
a decreases), the static Peierls stress increases and the dislocation becomes harder to move. 
The size of the dislocation core is also related to the shape of g(x). As shown by Figures 
El to the dislocation core expands as a increases: for a < 0.26, the dislocation core is 
very narrow, as shown in Figs. E^a) and Ufa). Figs. Efb), and Ufa) to c) show that the 
dislocation core widens one lattice point as a sweep the interval 0.26 < a < 0.29, in which 
the variation of the Peierls stress is very small (see the plateau in Fig. |2J). The dislocation 
core gains more lattice points as a increases to 0.32 and beyond, cf. Figs. Efc), b) and 
[Tfd) . Thus the size of the dislocation core and the Peierls stress are related to the width of 
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the interval of x for which g'{x) < and to the actual value of the slope. As an additional 
example, the symmetric sine function g(x) = sin(27nr)/(27r) has wide subintervals of small 
slope, which produces very low Peierls stresses and an artificially wide dislocation core. In 
|13l Il5| . this problem was avoided by setting w 2 = and changing g(Dfui) <-> D^Ui and 
27ig(D2Ui) = sin(2irDfui) in and (|18|) . A moving dislocation is a discrete traveling 
wave advancing along the x axis, and having far field (ui(l — ct,m) + Fm,U2{l — ct,m)). 
The analysis of depinning and motion of planar edge dislocations follows that explained in 
Ref. [3] with technical complications due to our more complex discrete model. 



C. Interaction of edge dislocations and crack formation 

The reduced 2D model (fTBj) - (j!7|) can be solved numerically to illustrate interaction of 
edge dislocations. Figure |S] illustrates the repulsion of equal-sign edge dislocations whereas 
opposite edge dislocations attract each other and form dislocation loops as in Fig. El or 
dislocation dipoles as in Fig. Friction terms affect numerical simulations of the model 
as follows. As in the case of ID models 18], atoms may oscillate far from the core of a 
moving dislocation when the equations of motion are conservative or slightly damped. Large 
friction (order one coefficients) reduces the oscillations of individual atoms, the instantaneous 
position of the core of the defect is easier to locate and its movement in the distorted lattice is 
easier to follow. Small friction (order 10 -2 coefficients) results in dislocation glide combined 



with oscillations o 
friction. See Refs. 



the individual atoms. The figures |H1 El and [TU] were obtained with small 



and 



19l on the impact of friction and inertia on ID wave front profiles. 
Figure ITT1 shows the formation of a crack propagating in the x direction under an applied 
tension in the y direction. In principle, numerically solving the discrete equations of motion 
we can find: (i) the threshold stress for crack propagation, (ii) the direction of propagation, 
(iii) the crack speed, and (iv) the crack shape. We do not need to impose additional con- 
ditions such as displacement thresholds for breaking atomic bonds as in the usual spring 



models for brittle fracture 
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FIG. 8: (Color online) Repulsion of like sign edge dislocations. 
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FIG. 9: (Color online) Attraction of opposite sign edge dislocations leading to formation of a 
dislocation loop. 

III. ELASTICITY IN A NON-ORTHOGONAL BASIS 
A. Equations of motion 

For fee or bee crystals, the primitive vectors of the unit cell are not orthogonal. To find 
a discrete model for these crystals, we should start by writing the strain energy density in a 
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FIG. 10: (Color online) Atraction of opposite sign edge dislocations leading to formation of a 
dislocation dipole. 
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FIG. 11: (Color online) Snapshots showing crack generation and growth induced by applying a 
tension in the direction y to a dislocation dipole. 

non-orthogonal vector basis, a±, 03, in terms of the usual orthonormal vector basis e\, e 2 , 
e% determined by the cube sides of length a. Let Xi denote coordinates in the basis e,, and let 
x\ denote coordinates in the basis dj. Notice that the Xi have dimensions of length while the 
x\ are dimensionless. The matrix T = (ai, a 2 , a 3 ) whose columns are the coordinates of the 
new basis vectors in terms of the old orthonormal basis can be used to change coordinates 
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as follows: 

%i ^ij TijXj. (20) 

Similarly, the displacement vectors in both basis are related by 

u 'i = T ij lu ji u i = T ijU'j, (21) 

and partial derivatives obey 

A = r ,A JL = T - l JL. ( 22) 

dx\ 3l dxj' dxi jl dx'j 
By using these equations, the strain energy density W = (l/2)ciki m eikeim can be written as 

2 ijlm dx j dx m 2 Crsm dx' s dx' q ' 1 } 

where the new elastic constants are: 

C rspq = c ijlmTi r T s j Ti p T qm . (24) 

Notice that the elastic constants have the same dimensions in both the orthogonal and the 
non-orthogonal basis. To obtain a discrete model, we shall consider that the dimensionless 
displacement vector depends on dimensionless coordinates x\ that are integer numbers 
= «•(/, m, n; t). As in the case of sc crystals, we replace the distortion tensor (gradient 
of the displacement vector in the non-orthogonal basis) by a periodic function of the cor- 
responding forward difference, tu,- = g{D^u' i ), cf. Eq. ©• The discretized strain energy 
density is 

W(l, m, n,t) = l - J rspq g{Dtu' r ) g(D+u' p ). (25) 

The elastic constants c' rspq can be calculated in terms of the Voigt stiffness constants for a 
cubic crystal, C u , C u and C 12 . Eq. © yields Cy 7m = C u (5 u 5 jm + 5 im 5ij) + C 12 5ij5 lm - 
H(8 li 5ij5 1 i5 lm + hi^2jhih m + S 3i 5 3j 5 3 i5 3m ), where H = 2C 44 + C 12 - C n measures the 
anisotropy of the crystal and Eq. (|2"4"j) provides the tensor c' rspq . The elastic energy can be 
obtained from Eq. (|25j) for W by means of Eqs. (0J). Then the equations of motion (JZj) are 

,3 



at 2 iq pq du' p 

which, together with Eqs. (JU) and (|23j). yield 

d 2 u' 

p-^ = Tr q ^Dj[c' pjrs g f (Dfu' p )g(Dtu f r )}. (26) 

This equation becomes (fTUj) for orthogonal coordinates, = 5i q /a, once we take into 
account the Einstein convention on summation over repeated indices in (JlOp . 
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B. Far field of a dislocation 

As in the case of sc crystals studied in Section [TT1 we should determine the elastic far field 
of a dislocation under zero stress to set up the initial and boundary data needed to solve 
numerically the discrete equations of motion (|26|) . We can calculate the elastic far field of 
any straight dislocation following the method explained in Chapter 13 of Hirth and Lothe's 
book Firstly, we determine the elastic constants in an orthonormal coordinate system 
e'(, e'2, e'l with e 3 ' = — £ parallel to the dislocation line. The result is 

3 

C ijkl = C ijkl H ^""^ {Sj n Sj n SknSln 8i n 5j n 5 kn 5l n ). (27) 
n=l 

Here the rows of the orthogonal matrix S = (e'(, e 2 , e 3 ')* are the coordinates of the e"'s in 
the old orthonormal basis ei, e2, e%. In these new coordinates, the elastic displacement field 
(u'l, u'2, U3) depends only on x'[ and on x'^. The Burgers vector and the elastic displacement 
field satisfy b" = b'l = and u'{ —u'2 — Q for a pure screw dislocation in an infinite medium. 
For a pure edge dislocation, 63 = and u'l = 0. Secondly, the displacement vector (u'{, u'^, u'l) 
is calculated as follows: 

• Select three roots Pi,P2,P3 with positive imaginary part out of each pair of complex 
conjugate roots of the polynomial det[a ik (p)] = 0, a ik (p) = c'{ lkl +(c'- lk2 +c'[ 2kl )p+c'- 2k2 p 2 . 

• For each n = 1, 2, 3 find an eigenvector A k (n) associated to the zero eigenvalue for the 
matrix a ik (p n ). 

• Solve ReE 3 n=1 A k {n)D{n) = b'l, k = 1,2,3 and ReELi ELi«2H + 
c'l 2k2 p n )A k (n)D{n) — 0, i — 1, 2, 3 for the imaginary and real parts of D(1),D(2),D(3). 

• For k = 1, 2, 3, u" k = M~m ^Li A k (n)D(n) ln(4 + p n x 2 % 

Lastly, we can calculate the displacement vector u' k in the non-orthogonal basis from u k . 

C. Discrete models for fee metals 

For fee metals, the non-orthogonal vector basis comprising primitive vectors is 

O! = ~ (1, 1, 0), a 2 = % (0, 1, 1), a 3 = - (1, 0, 1). (28) 
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The equations of motion are (|26|) with the corresponding transformation matrix T = 
(ai,a 2 ,a 3 ). 

We shall now analyze the motion of dislocations in the case of gold. The initial and 
boundary data for the numerical simulations are constructed from the far fields of disloca- 
tions in anisotropic elasticity as explained in the previous Subsection. We have considered 
two straight dislocations: the perfect edge dislocation directed along £ = (— 1,1,— 2)/\/6 
(with a Burgers vector which is one of the translation vectors of the lattice, and therefore 

n 

glide of the dislocation leaves behind a perfect crystal [9j) and the pure screw dislocation 
along £ = (1, l,0)/\/2. For the perfect edge dislocation, we select: 

4 = (-1, -1, 0)/V2, 4 = (1, -1, -1)/V3, 4 = (1, -1, 2)/V6, (29) 

which are unit vectors parallel to the Burgers vector b, the normal to the glide plane n, and 
minus the tangent to the dislocation line — £, respectively. For the pure screw dislocation, 
we have: 

4 = (1, -1, 2)/V6, 4 = (-1, 1, 1)/V3, 4 = (-1, -1, 0)/y/2, (30) 

where 4 i s a um t vector normal to the glide plane and e^' is a unit vector parallel to the 
dislocation line and to the Burgers vector (but directed in the opposite sense). 

For gold, Cn = 186 GPa, C 44 = 42 GPa, C 12 = 157 GPa and H = 55 GPa. The lattice 
constant is a = 4.08 A and the density is p = 1.74 g/cm 3 . Figures IT21 and UBI show the perfect 
edge dislocation and the screw dislocation obtained as stationary solutions of model (|26|). 
Due to the boundary conditions we have chosen, their far fields match the corresponding 
elastic far fields of the dislocations (written in the non-orthogonal coordinates corresponding 
to the primitive cell vectors a±, a^). Dark and light colors are used to trace points placed 
in different planes in the original lattice. Note that the planes perpendicular to the Burgers 
vector in Fig. ^] have a two-fold stacking sequence 'dark-light-dark-light . . . ' The extra 
half-plane of the edge dislocation consists of two half planes (one dark and one light) in 
the dark-light-dark-light . . . sequence. Movement of this unit dislocation by glide retains 
continuity of the dark planes and the light planes across the glide plane, except at the 
dislocation core where the extra half planes terminate jjj. 
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FIG. 12: (Color online) Perfect edge dislocation in a gold lattice displaying a two-fold stacking 
sequence of planes. Lines locating the dislocation core are a guide for the eye. 

D. Discrete models for bcc metals 

The discrete model for bcc metals is similar to that for fee metals explained in the previous 
Subsection, but the non-orthogonal vector basis comprising primitive vectors is now 

a x = | (1, 1, 1), a 2 = ~ (-1, 1, 1), a 3 = ~ (1, -1, 1). (31) 

The equations of motion are (|26j) with the corresponding transformation matrix T = 
(ai, a 2 ,a 3 ). 

As in the previous Subsection, we calculate the elastic displacements of an edge dislocation 
and a screw dislocation in iron. For the edge dislocation we select: 

e'i = (1, 1, 1)/V3, 4 = (-1, 0, 1)/V2, 4 = (1, -2, 1)/VQ, (32) 

which are unit vectors in the directions of the Burgers vector b, the normal to the glide 
plane n and the dislocation line vector, respectively. For the pure screw dislocation: 

e'[ = (-1,0, 1)/V2, 4 = (-1, 2, -1)/V6, 4 = -(1, 1, 1)/V3 (33) 

where e^' is the normal to the glide plane and 4 a unr k vector parallel to the dislocation line 
and to the Burgers vector. 
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FIG. 13: (Color online) Screw dislocation in a gold lattice. Lines locating the dislocation core are 
a guide for the eye. 

For iron, C u = 242 GPa, C 44 = 112 GPa, C 12 = 146.5 GPa and H = 129 GPa. The 
lattice constant is a = 2.87 A and the density p = 7.86 g/cm 3 . Figures IT41 and ITB1 show 
the edge and the screw dislocations obtained as stationary solutions of model (|26|). Their 
far fields match the corresponding elastic far fields of the dislocations (written in the non- 
orthogonal coordinates corresponding to the primitive cell vectors a\, a,2, 03). Dark and 
ligh colors are used to trace points placed initially in different planes perpendicular to the 
Burgers vector. 

IV. CONCLUSIONS 

We have proposed discrete models describing defects in crystal structures whose con- 
tinuum limit is the standard linear anisotropic elasticity. The main ingredients entering 
the models are the elastic stiffness constants of the material and a dimensionless periodic 
function that restores the translation invariance of the crystal and, together with the elastic 
constants, determines the Peierls stress. The parameter value of a specific one-parameter 
family of periodic functions can be selected so as to fit the observed or calculated value of 
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FIG. 14: (Color online) Edge dislocation in an iron lattice. Lines locating the dislocation core are 
a guide for the eye. 




FIG. 15: (Color online) Screw dislocation in an iron lattice. Lines locating the dislocation core are 
a guide for the eye. 

the Peierls stress for the material under study. For simple cubic crystals, their equations of 
motion are derived and solved numerically to describe simple screw and edge dislocations. 
Moreover, we have obtained numerically edge dislocation loops and dipoles, and observed 
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crack generation and growth by applying a tension in the vertical direction to a dislocation 
dipole. For fee and bee metals, the primitive vectors along which the crystal is transla- 
tionally invariant are not orthogonal. Similar discrete models and equations of motion are 
found by writing the strain energy density and the equations of motion in non-orthogonal 
coordinates. In these later cases, we have determined numerically stationary edge and screw 
dislocations. 
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APPENDIX A: DERIVATION OF THE EQUATIONS OF MOTION 



Firstly, let us notice that 

dW dW de jk 







dui(l,m,n;t) dej k dui(l,m } n;t) 



dui(l, m, n; t) 



d{D+u k ) 



g'(D+u k ) - , 1 K ', + g'(D+u 



g{Dfu k ) + g(D+ Uj 
d{Dtu 3 ) 



(Al) 



3 dui(l, m, n; t) a v fc 3 dui(l,m,n;t) 
where W is a function of the point (l',m',n'), and we have used the definition of stress 



tensor: 



dW 



(A2) 



and its symmetry, cr^- = a^. Now, we have 

d 

—— —[Dfu k (l', m, n'\ t)} = 5 ik {8u<+x - 8iv) S mm '5 r , 

oui{l,m,n;t) 

and similar expressions for j = 2, 3. By using (jAl|) - (jA3|) . we obtain 
d 

d Ui (l } m } n;t) v ±f 



(A3) 



^ W(V, m', n'; t) = - E D J[°H 9\D+u t )} . 



(A4) 



In this expression, no sum is intended over the subscript i, so that we have abandoned the 
Einstein convention and explicitly included a sum over j. Therefore Eq. (jHJ) for conservative 
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dynamics becomes 



Mui = J2 D 7K9'(Df Ui )}, (A5) 



3 



which yie.ds Ec, 03. Except for the factor ^>,), these equations are discrete versions 
of the usual ones in elasticity 16] . 
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